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ABSTRACT 

The stability characteristics of various compact fourth- and sixth-order spatial opera- 
tors are assessed using the theory of Gustafsson, Kreiss and Sundstrom (G-K-S) for the 
semi-discrete Initial Boundary Value Problem (IBVP). These results are then generalized to 
the fully discrete case using a recently developed theory of Kreiss. In all cases, favorable 
comparisons are obtained between G-K-S theory, eigenvalue determination, and numerical 
simulation. The conventional definition of stability is then sharpened to include only those 
spatial discretizations that are asymptotically stable (bounded, Left Half-Plane eigenval- 
ues). It is shown that many of the higher-order schemes which are G-K-S stable are not 
asymptotically stable. A series of compact fourth- and sixth-order schemes, which are both 
asymptotically and G-K-S stable for the scalar case, are then developed. 


‘This research was supported by the National Aeronautics and Space Administration under NASA Con- 
tract No. NAS1-18605 while the second and third authors were in residence at the Institute for Computer 
Applications in Science and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23665. 
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1. INTRODUCTION 


Recently, higher-order numerical methods have seen increasing use in the Direct Numeri- 
cal Simulation (DNS) of the Navier-Stokes equations. Although they do not have the spatial 
resolution of Spectral methods, they offer significant increases in accuracy over conventional 
second-order methods. They can be used on any smooth grid, and do not have an overly 
restrictive CFL dependence as compared with the 0(N~ 2 ) CFL dependence observed in 
some Spectral methods on finite domains. In addition, they are generally more robust and 
less costly than Spectral methods. The issue of the relative cost of higher-order schemes 
(accuracy weighted against physical and numerical cost) is a very complex issue, ultimately 
depending on what features of the solution are sought and how accurately they must be 
resolved. In any event, the further development of the underlying stability theory of these 
schemes is important. 

The state of higher-order temporal discretizations is well developed and relies greatly 
on the existing Ordinary Differential Equation (ODE) literature. Higher-order spatial dis- 
cretizations are well documented in the literature, as well. For example, entire classes of 
centered explicit spatial schemes are described in the text of Kopal [1]. In practice, compact 
schemes (methods where both the solution and its derivatives are treated as unknowns and 
are solved for simultaneously) are more accurate than optimal explicit schemes (nth order 
schemes involving n+1 grid-points), and have gained much attention recently for use in DNS. 
The fundamental ideas of compact schemes as well as derivation techniques can be found in 
the work of Vichnevetsky [2], and will not be pursued in this work. The primary difficulty 
in using higher-order schemes is finding stable boundary schemes that preserve their formal 
accuracy. It is well known that for a hyperbolic system to preserve formal spatial accuracy, 
an ( N) th order inner-scheme must be closed with at least an (N—l) th order boundary scheme 
[3]. To date, many higher-order inner-schemes are used with lower-order boundary schemes, 
because no stable higher-order formulations are known. The formal accuracy of these for- 
mulations is thus reduced to one order more than the Boundary Condition (BC) accuracy, 
and it is questionable whether the additional work incurred in the higher-order inner-scheme 
is justified. 

Determining the numerical stability of a fully discrete approximation for a linear hy- 
perbolic partial differential equation is a difficult task. For the “Cauchy Problem” on an 
infinite domain (— 00 , 00 ), standard techniques based on Fourier methods generally provide 
the necessary conditions for stability of the numerical scheme. For the Initial-Boundary- 
Value Problem (IBVP) on the semi-infinite domain [0, 00 ), or the finite domain [-1,1], Fourier 
techniques are not straight forward to apply and do not provide sufficient conditions for nu- 
merical stability. To address these issues, Osher [4] and Kreiss [5], and later Gustafsson et 
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al. [6], developed stability analysis techniques based on normal mode analysis. Their work 
(generally referred to as G-K-S stability theory) established conditions that the inner and 
boundary schemes must satisfy, to ensure stability. The G-K-S theory states that hyperbolic 
systems are assured stability if no eigenvalues or generalized eigenvalues exist for the IBVP. 
Trefethen [7] further clarified the physical meaning of the G-K-S condition by noting that 
the concept of stability at a boundary could be related to the group velocity of the boundary 
scheme, specifically, whether it is carrying energy into or out of the numerical domain. 

One of the weak points of fully discrete G-K-S theory has been the great complexity 
in applying it to higher-order numerical schemes. Raising the spatial or temporal accuracy 
generally increases the complexity of the stability polynomials (the order increases, giving 
rise to more roots) which govern the stability of the numerical approximation. In multi-stage 
time discretization schemes (e.g. Runge-Kutta schemes with three or more stages), where 
boundary conditions must be applied at intermediate levels, the stability polynomials that 
must be tested at each boundary are nearly insurmountable using analytic techniques. One 
can greatly simplify the analysis by addressing the semi-discrete problem as a method-of-lines 
IBVP, rather than the fully discrete problem. The underlying G-K-S theory for the semi- 
discrete problem was initially developed by Strikwerda [8], He showed that by discretizing 
space and leaving time continuous, the necessary and sufficient conditions for method-of-lines 
IBVP stability are analogous to those governing the stability of the fully discrete case. 

The precise connection between the semi-discrete stability bounds and those obtained in 
the fully discrete analysis is not always straight forward. Recently, however, Kreiss et al. [9] 
have shown that under very weak conditions, stability of the semi-discrete approximation 
infers stability in the fully discrete approximation if specific Runge-Kutta time marching 
schemes are used. Therefore, one can rely on the semi-discrete G-K-S theory to assess the 
stability with R-K integration, of various higher-order spatial discretization operators, thus 
simplifying the calculations considerably. 

The emphasis of this work is to apply semi-discrete G-K-S theory to several higher- 
order spatial operators for the IBVP. Since compact methods naturally lend themselves 
to fewer implementational difficulties at the boundaries, they will be the primary focus. 
Stable boundary formulations which preserve the formal accuracy of the inner-scheme will 
be presented for spatial derivative operators of up to sixth-order. 

2. MODEL EQUATION FOR IBVP 

Under the assumption of locally frozen coefficients, the equations governing conservation 
of mass, momentum and energy in one spatial dimension can be transformed into a system 
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of hyperbolic equations, having the form: 


dU . dU 

~dt ~ A lh +F] x ~ °* * - 0 


(i) 


where “A” is a diagonal matrix with real eigenvalues, F is a source term, and boundary and 
initial data are supplied in the form 


t/ 7 (0, t) + TU n (0, t ) = g(t), U(x, 0) = f(x) 


( 2 ) 


where T is a matrix describing the boundary conditions, and A has been divided into its 
right- and left-going characteristics to determine the boundary conditions. The problem is 
said to be well-posed if the solution U(x,t) depends smoothly on the initial and boundary 
data. Our focus in this work will be on the scalar form of eqn (1-2), where the matrix “A” 
is a negative real constant “a” and with the source term F = 0. This simplification can 
be justified since stability of a numerical scheme on the scalar equation implies stability of 
the system if boundary conditions are imposed in a characteristic form (ex. see Gottlieb et 
al. [10]). A further necessary modification is that the semi-infinite spatial domain must be 
truncated to a finite domain. For simplicity, we shall assume that the physical domain is 
confined to the interval 0 < x < 1 for all times. 


3. SEMI-DISCRETIZATION 

In this work, the numerical discretization of eqn (1) will be accomplished by two separate 
and independent steps. The spatial derivatives will first be approximated with appropriate 
formulas, leaving what is generally referred to as a semi-discretization. The numerical solu- 
tion will then be advanced in time in a Method-of- Lines approach, using a stable temporal 
scheme. We begin by dividing the continuous domain [0,1] into N uniform intervals of width 
Ax where N Ax = 1 . The continuous derivative U x is then replaced with a finite-difference 
representation involving the functional values Uj at the discrete points. A system of ODE’s 
results having the form 

dV, 

dt 

where 


= M + Vj j = 0, ... , N 


( 3 ) 


M + Vj = £ a m*V j+k 

k=-L, 


( 4 ) 


and “L” and “R” are the width of the stencil extending to the left and right of grid-point 
j, respectively. Note that m^Lj and Rj are functions of “j” since there is no reason to 
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assume that the same stencil will be used at each grid-point. For example, consider the 
case where the same spatial stencil is used at every interior grid-point in the domain. We, 
therefore, choose a particular Lj = L and Rj = R. (Obviously, the choice of L, R, and 
m k - should coincide with a scheme that is stable for the Cauchy problem on the infinite 
domain). This scheme can only be used for L < j < N — R without the stencil protruding 
through the boundary. Thus, exactly L + R additional formulas have to be defined near the 
boundaries. Since there is only one physical boundary condition, L + R - 1 of the schemes 
are strictly numerically motivated. These schemes are generally referred to as numerical 
boundary schemes (NBS). 

Noting that the physical boundary condition g(t) will be imposed at the grid-point j=0, 
eqn (3) can be rewritten as 

dV 

= MVj + Bjg(t)-, j = 1 N (5) 

where M is an N x N matrix, and B } is a vector of dimension N, describing the dependence 
of the “j th ” scheme on the boundary data. The matrix M usually is diagonal of order L + R 
+ 1 for most explicit methods, but can in general be full. With little loss of generality, g(t) 
can be chosen to be zero. The exact solution of the semi-discretization described by eqn (5) 
for homogeneous boundary data becomes 

V 3 (t) = /(xj)ex p Mt ; j = 1, ... ,N. (6) 


Note that all the boundary information is incorporated directly into the matrix M. and that 
the stability of the numerical scheme depends directly on the properties of the matrix, not 
just on what spatial discretization was chosen for the inner-scheme. 

It is instructive to clarify these points involving boundary closures with examples of an 
explicit and an implicit spatial scheme. We choose to concentrate on schemes having at 
least fourth-order spatial accuracy since most of the difficulties associated with high-order 
stencils are not observed with second-order schemes. The first is an explicit 5-point scheme 
reported by Gary [11], and later shown by Strikwerda [8] to be G-K-S stable. The scheme is 
uniformly fourth-order in space. The spatial discretization is accomplished with the stencils 


dVo 

dx 

dVi 

dx 

dx 

dV N - i 

dx 

dV n 

dx 


1 

12Ax 

1 

12Az 

1 

12Ax 

1 

12Ax 

1 

nKx 


(—25 Vo + 48Vi - 36V 2 + 16V 3 - 3V 4 ) 

(— 3Vo - 10Va + 18V 2 - 6V 3 + V 4 ) 

(Vji — 2 - 8V,-! + SVj+i - v j+2 y, J = 2, 

(— V)v_4 + 6 Vat-3 — I 8 V/V -2 + 10 Vyv-i + 3VV) 

(3V/v_4 - l6V N - 3 + 36Vyv_ 2 - 48VJv-i + 25V N ). 


N — 2 


( 7 ) 
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The terms in eqn (3) take the form 


Vo ' 


' -25 

48 

-36 

16 

-3 







V 


-3 

-10 

18 

-6 

1 







V 2 


1 

-8 

0 

8 

-1 




0 



V 3 



1 

-8 

0 

8 

-1 






■ 

, M+ = 



* 


. 

* 

• 

. 



(8) 

V N -2 







1 

-8 

0 

8 

-1 


V N . i 


i 


0 



-1 

6 

-18 

10 

3 


V N 







3 

-16 

36 

-48 

25 



where Note that four NBS’s are required to close the numerical scheme, of 

which only one will be overwritten by a physical boundary condition. Accounting for the 
boundary condition at the grid-point j=0 in eqn (8), results in expressions for eqn (5) of the 
form 


' V, ' 


’ -10 18 -6 1 

V 2 


-80 8-1 

V 3 


1-808 

■ 

, M — 

• 

Vn-2 



Vjv-i 


0 

. V N 


L 




' -3 ' 

0 


1 

-l 


0 

. 

, B = 

• 

1 -8 0 8 -1 


0 

-1 6 -18 10 3 


0 

3 -16 36 -48 25 


1 

o 


( 9 ) 


where M = and B = j^B. It is evident that the matrix M is the N x N submatrix 

of the matrix M + , corresponding to 1 < i,j < N. 


The implicit fourth-order example is that of a compact discretization in which the nu- 
merical approximation to ^ is made in the following form 


dVo JW 

dx + dx 

av,-, dVj av i+ , 

dx dx dx 

dV N -i dV N 

dx dx 


^j(-17K„ + 9V, + m - v 3 ) 
g^(Viv-3 — 9V1 v- 2 — 9Vyv-i + 17 Vat). 


( 10 ) 


Noting the structure of the spatial operator and the fact that ^ for this scheme, 

eqn (3) is often rearranged into the form 


dVi 


M '-df = M ? Vii 3= 0, ...,N 


( 11 ) 
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where M + = (M x + ) 1 M 2 . The resulting matrix expression becomes 



i 

^ ST 


' 6 18 
1 4 1 

0 

y+ = 

• ■ 

l 

, Mt = 

0 

1 4 1 

18 6 _ 


Mo = 


-17 9 9 -1 

-3 0 3 0 


0 -3 0 3 

1 -9 -9 17 


( 12 ) 


where M 2 = Imposing the boundary information at grid-point j=0 reduces the 

matrix by one order and yields matrix equations of the form 


M 2 - 



V\ 



6 6 




v 2 



1 4 1 


0 

V = 

• 

5 

Mi = 

• 

■ • 



V N -! 



0 

1 

4 1 


Vjv 



- 


18 6 _ 

'-19' 

9 




' -1 ' 


-3 0 : 

3 

0 



0 






b 2 = 

• 


0 

| 

3 0 3 



0 



1 -' 

9 -9 17 _ 



0 



(13) 


where B 2 = ^f? 2 . Note that because a matrix multiplication was involved in determining 
the spatial operator for grid-points 1 < j < N, matrix M is not a simple submatrix of matrix 
M + as it was with the explicit scheme. The spatial operator described by eqn (13) is referred 
to as being compact because both M x and M 2 only involve three grid-points (except near 
the boundary in M 2 ). Because of this compactness, L = R = 1 , and only two NBS must 
be defined (one of which includes, but is not replaced by, the physical boundary condition). 
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Expressing this algorithm in the form of eqn (5) yields 

M = M^M 2 , B = M~ X B 2 . (14) 

It is apparent that while the algorithm is implementationally compact, the resulting M is a 
full N x N matrix. 

Working with the scalar form of eqn (1), with the boundary conditions posed at x=0, 
has allowed us to simplify the semi-discretization defined by eqn (3) to that of eqn (5). Note 
that if the governing equation would have been Ut + aU x = F with boundary data at x 
— 1, the two previously mentioned schemes could have easily accommodated this change 
and still produced a stable algorithm. For a system of hyperbolic equations with variable 
coefficients, one does not know a priori the sign of the eigenvalues of the matrix A in eqn (1). 
Therefore, the solution is advanced in time for j=0,N, followed by a characteristic decom- 
position in which the physical boundary conditions are imposed at either j=0 or j=N such 
that the resulting system is well-posed. This requirement imposes a constraint on the class 
of allowable matrices M + which can be conveniently used to obtain stable discretizations. 
Only central-difference schemes have this desirable feature; that of being stable for a > 0, or 
a < 0. Assuming that a central-difference operator is used throughout the interior domain, 
then eqn (4) requires that R = L NBS’s be defined at each boundary. The structure of these 
NBS’s must be such that the resulting scheme be stable for either the inflow (a < 0; j = l,N) 
or the outflow (a > 0; j=0,N-l). This can only be accomplished if the same methodology is 
used to derive the NBS s at each end of the domain. Therefore, the NBS’s at each respective 
boundary, are asymmetric mirror images of each other. The resulting matrix M + has the 
following property 


M + = —PM + P (15) 

where P is the permutation operator defined by = 0,0 < i,j < N , for i ^ N - j, and 

Pi,j = TO ^ i,j < N , for i = N — j. Note that PP = I, the identity matrix. The matrix M 
is, therefore, the N x N submatrix of M + . 

We now have a well defined class of spatial operators which are acceptable for our pur- 
poses. They are high-order central difference schemes with boundary implementations which 
are imposed asymmetrically. For future reference, we shall define a convenient nomencla- 
ture to describe these matrices. The matrix M + shall be described as (NBS U - • -,NBS L - 
CD — NBSr, • • • NBSn ), where “CD” is the order of the central difference operator used 
in the interior, and NBSj is the order of the NBS used to close the scheme at the points 
next to the boundary. For example, the explicit uniformly fourth-order scheme represented 
by matrix eqn (8) is denoted by the nomenclature (4, 4-4-4, 4), where the “-4-” denotes the 
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inner-scheme being approximated with a fourth-order stencil, and the symmetric “4,4” de- 
notes fourth-order stencils at j=0,l, and j=N-l,N. The three-point compact scheme described 
by eqn (12) is denoted (4-4-4). Again, the “-4-” refers to the inner-scheme accuracy, and 
the symmetric “4”s indicate closure on the boundaries of fourth-order accuracy. There is 
no ambiguity in the nomenclature between the compact and explicit schemes, since only the 
compact scheme can retain fourth-order inner accuracy with one NBS at each boundary. 

4. STABILITY OF THE IBVP 

The eigenvalues of matrix M from eqn (5), resulting from higher-order finite-difference 
approximations to U x , tend to align along the imaginary axis in complex conjugate pairs. 
To time advance eqn (5) efficiently, the time discretization algorithm should include in its 
stability domain a large portion of the imaginary axis. Conventional Runge-Kutta (R-K) 
time advancement algorithms of third- or fourth-order are well suited for semi-discretizations 
of hyperbolic equations, and are the only method considered in this study. In particular, 
the standard fourth-order method of Kutta (ex. see Gear [12]) will be used in this study 
because of the 1) fourth-order non-linear accuracy, 2) the large stability envelope, and 3) the 
low storage requirements. 

It is desired to know whether a given higher-order spatial discretization is stable for this 
time advancement scheme. Using conventional G-K-S theory for the fully discrete IBVP for 
fourth-order R-K time, and second-order space discretization would involve polynomials of 
eighth-order in k to solve at the boundaries. Closed form solutions would be difficult to 
obtain under these circumstances, and would get more complicated with increasing spatial 
accuracy. The stability analysis can be greatly simplified by relying on three fundamental 
theorems of stability analysis which are valid under the conditions proposed in this study. 
We will discuss each briefly, but for further clarification suggest consulting the original works. 
The essential elements of the theorems forming the basis of this work are: 

Theorem 1: G-K-S theory (fully discrete [6] or semi-discrete [8]) asserts that to show 
stability for the finite domain problem, it is sufficient to show that the inner-scheme is 
Cauchy stable on (-oo, oo) , and that each of the two quarter-plane problems is stable using 
normal mode analysis. Thus, the stability of the finite-domain problem is broken into the 
summation of three simpler problems. 

Theorem 2: For each quarter plane problem arising in theorem 1, a necessary and 
sufficient condition for stability of the initial boundary value problem is that there exist no 
eigensolution. This theorem is true for either the fully-discrete case [6], or the semi-discrete 
case [8]. 
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The algebraic complexity involved in showing stability of the IBVP is dramatically re- 
duced in the semi-discrete case, since time remains continuous. Ultimately, numerical stabil- 
ity is a fully-discrete concept, and a connection between the semi- and fully-discrete stability 
must be used. The third theorem provides this connection. 

Theorem 3: Under mild restrictions [9], if a semi-discrete approximation is stable in a 
generalized sense and a Runge-Kutta method which is locally stable is used to time march 
the semi-discretization, the resulting totally discrete approximation is stable in the same 
sense so long as the stability region of the R-K method encompasses the norm of the semi- 
discretization. 

It should be noted [9] that the stability definitions used in the first two theorems (G-K-S 
stability) are different than that of the third (generalized stability). The first two theorems 
rely on G-K-S stability (sometimes referred to as the Kreiss condition) in the semi- or fully- 
discrete case. At least two different definitions of G-K-S stability are encountered in the 
cumulative works on G-K-S analysis. They are: 

Definition 1: The IBVP is stable if for r) > rj 0 , the solutions of eqn (1) with homogeneous 
initial data satisfy 

|i>(0, S)f + (, - ,„) 2 || U(., S)|| 2 < K(\ g\ l + ||F(., S)|| 2 ) (16) 

where r) 0 and K are universal constants, and the "denotes the Laplace transformed variables. 

Definition 2: The IBVP is stable if for r] > r] 0 , the solutions of eqn (1) with homogeneous 
initial data and g = 0 satisfy 

( n -r, 0 )\\U(.,S)\\<K(\\F(.,S)\\). (17) 

Both of these definitions can be related to what is often referred to as Lax stability, which 
is the numerical equivalent to an energy estimate satisfied by the discrete form of eqn (1,2). 

Theorem 3 relies on a different definition of stability that can be expressed as 

Definition 3: The IBVP is stable if for homogeneous boundary data (g=0) an estimate 
of the form 

||(/(.,f)|| = /Cexp“<‘-*>(||(7(.,t e )|j + /' ||F(.,t)||*-). (18) 

Jo 

Definition 2 is the most restrictive of the three conditions, but it has been conjectured 
by Kreiss [9] that Definition 2 and 3 are equivalent. The subtle differences between these 
definitions of stability will not affect the conclusions in this work, and the terms G-K-S 
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stability and Lax-stability shall be used interchangeably. We now describe in some detail 
the implications of the three Theorems. 

The chief difficulty in stability analysis is not which definition is chosen, but rather 
applying that definition to practical numerical methods and obtaining a stability bound. 
Since Fourier methods are not well suited on the finite domain, stability analysis can be 
carried out by energy methods or by normal mode analysis. Energy methods are, in general, 
very difficult to perform on high-order schemes. The modal relationships are simple to define 
but analytic solutions are often intractable. 

Theorem 1 describes how G-K-S analysis can be used to augment finite domain modal 
analysis. The original finite domain modal analysis is broken into the analysis of three 
equivalent, yet simpler modal problems. Assuming that a Cauchy stable scheme is used for 
the interior grid-points, the inner-scheme is tested for stability at each boundary in a semi- 
infinite spatial domain. In so doing, the stability of each boundary is independent of the 
influence from the other boundary. Stability of the two boundary problems implies stability 
of the finite domain numerical method. In addition, Theorem 1 provides a “perturbation 
test” for “generalized eigensolutions.” The test establishes the stability of certain borderline 
cases in normal mode analysis. 

To fully appreciate the power of G-K-S analysis, a normal mode analysis of the fourth- 
order compact scheme (4-4-4) described by eqn (13) is presented for the coupled finite domain 
problem. We proceed by assuming that the semi-discrete problem defined in eqn (11) has a 
solution of the form 


Vj(t) = exp 5 * <f>j (19) 

where “S” are the eigenvalues of the matrix M+ 1 . Substitution into eqn (11) yields the 

generalized eigenvalue problem 


MtSV 3 = M ? V (20) 

(for which we have assumed g(t) = 0), and the resolvent equation provided by the inner- 
scheme is 


{<f>j - 1 + + <f>j+\)S — 3(— <f>j-\ + <f>j+ 1) j — 2, ... , N — 1 (21) 

where S = It is apparent that <f>j = (j) 0 K j will satisfy this expression, yielding the 

following equation for the eigenvalues 

(t + 4 + K)5 = 3(— + k). (22) 

K K 
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This is a quadratic expression in k, and there will, in general, be two solutions which will 
satisfy eqn (22). The eigenvectors are, therefore, of the form (f>j = C\K + C 2 k 2 ^ . Simple 
manipulations show the two roots are related by 


and 


K i = K, K 2 


— 2 — K 

2 k T 1 


(j)j — C\ + C 2 ( 


-2-/c j 
2 k -(- 1 


We note that if k = a + ib, then 


(23) 


(24) 


1*1 1 — O? + fe 2 , 


*2! = 


b 2 + (a + 2) 2 
(26) 2 + (2 a -f 2) 2 


(25) 


and \ Kl \ > 1 implies |k 2 | < 1. Thus, each k is dominant near one of the boundaries and its 
influence decays monotonically with increasing distance from that boundary. The form of C\ 
and C 2 must be determined from the boundary conditions at j=l and j=N. The conditions 
at the two boundaries can be represented as 


C'i[S'-F 1 (kx)] + C 2 [5-F 1 (k2)] = 0 

Ci[§ - F n (k x)] 4- C 2 [S - F^(k 2 )] = 0 (26) 


where F\ and F )v are the functional relations resulting from substituting the eigenvector 
and eigenvalue into the expressions at grid-points j = l and j=N, respectively. Note that 
S = S(k), and that every term is a function of k. It is apparent that no general solution to 
this problem exists except the trivial one. Non-zero solutions exist for the condition where 
the determinant is equal to zero. The determinant condition gives the resulting expression 
for «, the roots of which, with eqn (22), give the eigenvalues and eigenvectors of the system. 

No closed form expressions are known for the roots to the determinant polynomial in 
this case or for any other reasonable boundary closures. The numerical solution of the 
determinant polynomial involves finding the roots of an “n th ” order polynomial in k. It 
is apparent that the boundary conditions dramatically influence the eigenvalue spectrum 
in the matrix M. Only in the limit N — ► 00 do the boundary conditions decouple. The 
power of G-K-S theory comes from breaking the normal mode problem into three separate 
problems. The roots to the k polynomials do not depend on the boundary to boundary 
coupling prevalent in eqns (25-26). 

Theorem 2 describes what constitutes stability for the two IBVP’s in Theorem 1 for the 
fully- or the semi-discrete case, and states that eqn(l) must satisfy the condition that no 
eigenvalues or generalized eigenvalues should exist for 7?.(S) > T) 0 . 
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Both theorems 1 and 2 rely on a definition of an eigensolution for their quarter-plane 
analysis. Here, an eigensolution is presented for the semi-discrete case. Similar definitions 
exist for the fully discrete case. Referring to the semi-discretization defined by eqn (5), 

Definition 4'- An eigensolution for the IBVP defined by eqn (5) is the nontrivial function 
V(x,s) satisfying [8] : 

I SV = M V x > 0 

II V J (0,s) - TV /7 (0,s) = 0 

III K(S) > 0 

IV for 'R-iS) >0 , V(x,s) is bounded as x — > oo 

V for 7Z =0 and |k| = 1 , a perturbation inside the unit circle of « (|k| = 1 — e, e >0) 
cannot produce an eigenvalue 7 Z> 6,8 >0. 

Note that since V’ = V’o*' 7 ? condition (IV) implies |/c| < 1. 

We refer to an eigensolution of the form (IV) or (V) as a G-K-S eigenvalue or a generalized 
G-K-S eigenvalue, respectively. With these conditions, the test for numerical stability has 
been simplified from the coupled normal mode analysis to tests involving 7l(S) > 0, for 
|/c| < 1 at each boundary, plus the exceptional case when both H{S) = 0, and |/c| = 
1. Theorem 3 relates the stability of the semi-discretization to the stability of the fully 
discrete numerical method. Obviously, analogous stability definitions must be chosen for 
both the semi- and fully-discrete cases, specifically the generalized stability condition given 
by Definition 3. Theorem 3 relies on temporal advancement schemes which are locally stable 
numerical methods. For a locally stable numerical method, the stability domain \z\ < 1 
(z is the amplification factor) in the complex plane, encompasses within the LH-P an open 
semi-circle of radius “Rj” centered at the origin, and symmetric about the real axis. The 
standard fourth-order R-K method satisfies this condition. Discretizing time with the fourth- 
order R-K method in eqn (5) produces a fully discrete method defined by 

V j {t + k) = L{kM)V j {t) + L{k)B J g{ty, j = 1, ...,iV (27) 

where the time-step A t = k , and where u L(kM)” is the polynomial in “k M” describing 
the time discretization. Then under very mild restrictions (see Kreiss [9]) on the eigenvalue 
structure of the matrix “L(kM)”, if the semi-discrete approximation is stable in a generalized 
sense, the totally discrete approximation is stable in the same sense for ||A:M|| < R\. 
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We have outlined a systematic approach to addressing the finite domain stability problem 
for the fully-discrete numerical approximation to eqn (1). The remainder of this work will 
describe the application of these techniques to several higher-order finite difference schemes. 


5. FOURTH-ORDER BOUNDARY CONDITIONS 


Before proceeding with the stability analysis of various higher-order boundary conditions, 
an example is presented that illustrates the necessity of (N - l) th order boundary closure 
for an N th order inner-scheme. The example also provides a numerical test to verify that 
the G-K-S theory is accurately predicting the stability behavior of the various numerical 
schemes. Consider the method-of-lines approximation to the scalar wave equation 


dU dU 

rH ~ O’ — 1 ^ x < 1) t > 0 


(28) 


U{t, — 1) — sin27r(— 1 — t), U(0, x) = sin27ra: — 1 < x < 1, t>0 (29) 

where the spatial discretization is accomplished by the fourth-order compact scheme de- 
scribed in detail by eqn (12). The exact solution is 


U(t,x ) - sin27r(a; — t), — 1 < x < 1, t > 0. 


(30) 


The comparison of the exact solutions with that obtained from various lower-order closures 
to the fourth-order inner-scheme provides a measure of the boundary influence. We begin 
by showing a grid convergence study performed to show the formal accuracy of the resulting 
schemes. The boundary condition formulas expressed at the grid-point “j=0” were 


9V 0 dVi —5V 0 + 4Vi + V 2 

~d^ + 2 ~d^ = 2Ax < 31 ) 

dVo^dVx _ -Vo + K 
dx dx 2Ax 

dV 0 -V 0 + Vt 

dx Ax ( 33 ) 

which represent third-, second-, and first-order closure at the inflow boundary, respectively. 
As mentioned earlier, the physical boundary condition was imposed at the point “j=0,” and 
the actual closure occurs at the point “j = l”. The closure could have been written explicitly 
for the point “j=l” by combining these formulas with the inner-scheme at “j = l” 


Wo A dV_ m — 3Vp + 3U 2 

dx dx dx Ax 


(34) 
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resulting in formula of the form 


dVi dV 2 

2 dx + dx 
dV, dV 2 

dVi dV 2 

4 — - H 

dx dx 


-Vo- 4Vi + 5V 2 

2Ax 

- Vo - 2Vi + 3V 2 

A T 

— 2V / o — Vi -4- 3 V 2 

Ax 


At the outflow boundary, closure was accomplished with the expressions 


9Vn dVN-i 
dx dx 

dVN gVfr-i 

dx dx 

dV N 

dx 


— 5Vv + 4Vjy— x + Vn- 2 

2Ax 

— Vn + Vv-i 

2Ax 

— Vn + Vjv-i 

Ax 


(35) 

(36) 

(37) 

(38) 

(39) 

(40) 


which represent third-, second-, and first-order spatial accuracy, respectively. These expres- 
sions are valid for the point “j=N” since there is no physical condition to impose there. 

In all cases, the temporal discretization was accomplished with the fourth-order Runge- 
Kutta algorithm. At every iteration, the solution was advanced for the grid-points j=0,N , 
followed by the imposition of the boundary condition at the point j=0. Since the inflow 
boundary condition was a nonlinear function of time, care was taken to evaluate the function 
corresponding to the proper intermediate level in time. Failure to do so degraded the formal 
accuracy of the method. The CFL used in the simulations was in the range 0.1 < CFL < 1. 
In no case did this violate the Von Neumann stability condition for the Cauchy problem. It 
should be noted that the formal truncation of the method is 0{At 4 ), and the error in time 
decays to the fourth power of the CFL for a given grid, and can be made as small as desired 
by decreasing the CFL. Typically, CFL’s < 0.1 were not needed since the dominant terms 
in the modified equation were negligible compared with spatial terms. Further reduction of 
the CFL resulted in no change in the error of the scheme. Recognition of this fact allows one 
to test the formal accuracy of methods with spatial accuracy higher than fourth-order for 
sufficiently small values of the CFL, and was used in the sixth-order simulations to determine 
formal accuracy. Finally, a third-order Runge-Kutta was used to test the generality of the 
temporal discretization in several cases. The results were quantitatively similar. 

The simulations were all run to equivalent times T=25, for all grids and methods at a 
CFL of 0.25. The error at T was then calculated and reported as an L 2 norm. The L norm 
produced similar results, but is not reported here. For methods which are Lax-stable, the 
error is bounded uniformly at each stage for 0 < A t <t and 0 < kAt < T, where k is the 
number of time-steps. Doubling the grid at constant CFL, should decrease the error at time 
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level T by a factor 2 where p is the order of the method. The formal accuracy of each 
scheme was determined for each of the Lax-stable schemes in this manner. 

Figures 1-3 show the results of this study. The log 10 of the L 2 error is plotted as a 
function of the log 10 of the number of grid-points resolving one period of the sine wave. 
The grid density ranges from 10 to 25 grid-points/(27r radians). Note that once a threshold 
accuracy is attained the points appear to decrease linearly. The slope of the data (dY/dX) 
gives the apparent formal accuracy of the method. In all cases, formal fourth-order accuracy 
is obtained with third- or fourth-order boundary conditions. Second order closure results 
in third-order formal accuracy, and first-order closure results in second-order accuracy. It 
is also apparent that the accuracy is greater for the fourth-order than for the third-order 
closure, regardless of the formal accuracy. These results are in agreement with the theory of 
Gustafsson [3] predicting that to retain formal accuracy of a numerical method, boundary 
conditions of order (N-l) must be imposed to retain N th order. Another interesting feature is 
that imposition of lower-order boundary conditions at the outflow plane results in a greater 
degree of error than at the inflow plane. One might have thought that since the characteristic 
is pointing out of the domain at the outflow boundary, error would be swept immediately 
out of the domain. This is obviously not the case. 

We will now derive the formal stability of the numerical boundary conditions used in 
this example. It can be shown that the fourth-order compact scheme is “Cauchy stable” 
for CFL < 1.63. I he stability of the inflow and outflow boundary conditions on the semi- 
infinite domain must be demonstrated. We begin by testing the outflow stability. The partial 
differential equation is 


dU dU 

_.___ = 0 , z>0, <>0. 


(41) 


No boundary condition is required in this problem, although a NBS is imposed at x = 0. As 
was done on the normal mode analysis for the finite domain, we assume a solution of the 

form Vj{t) = exp 5t <t> } , where <f>j = <f> 0 n j . Substitution into the inner scheme produces the 
resolvent condition for the eigenvalue S 


( — h 4 -f- k)S --- 3( f- /c) 

K /C 


(42) 


at each grid-point j > 1. At grid-point j=0, the scheme was closed with one of the bound- 
ary expressions given in eqns (12, 38, 39, 40) (obviously written for the grid-point j=0). 
Substitution of Vj(t) into these expressions yields 


(6 + 18 k)S = 
(2 + 4k)5 = 


— 17 + 9ac + 9k 2 — k 3 
— 5 + 4k + k 2 


(43) 

(44) 


15 



= 0 

(4 t/l ) 

(47) 

= 0 

(3 rd ) 

(48) 

= 0 

(2 nd ) 

(49) 

= 0 

(1 st )- 

(50) 


(1 + k)S = -2 + 2k (45) 

S = -1 + k. (46) 

Solving eqn (42) for S and substituting into eqns (43, 44, 45, 46) yield polynomials in k of 
the form 

(«-l) 5 

(k 2 + 4k + 1) 

(k 2 + 4k + 1) 

(*~1) 3 
(k 2 + 4k + 1) 

(k- 1) 2 (k + 2) 

( K 2 + 4/C + 1) 

Clearly, the only value of k which will simultaneously satisfy the inner resolvent condition 
and any of the outflow boundary conditions is k — 1. Substitutions of k = 1 into the 
resolvent condition produces 5 = 0, the case for which the perturbation test in G-K-S 
theory (condition V) must be used to show stability. Substituting k = 1 — e into the 
resolvent condition produces to first order 

65 = 6c. (51) 

For e > 0 and |/c| < 1, 5 < 0 showing stability of the perturbation. All the tested outflow 
boundary conditions are G-K-S stable, dhus, k = 1 is not a generalized eigenvalue. 

To show stability of the inflow conditions we study the partial differential equation 

^ ^ = o, X > 0, t > 0; (52) 

ot ox 

with the boundary condition imposed at x = 0 of the form 

U(0,t)=g(t), t> 0. (53) 

In spite of the physical boundary condition being imposed at j=0, a NBS must be imposed 
at j=l, and must be tested for stability. Substituting V 3 (t) = exp St <j>j, where <j>j = , into 

the inner-scheme produces the resolvent condition identical to eqn (42). Substitution into 
the first through fourth-order boundary conditions defined by eqns (12, 35, 36, 37), produce 
equations for S of the form 


(6 + 6k)5 = 

— 9 -f- 9 k -)- k^ 

(54) 

(4 + 2 k)S = 

—4 + 5k 

(55) 

(6 + k)S 

—4 + 6k 

(56) 

i 1 + 2 k)S — 

—2 + 6k. 

(57) 
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Note that without loss of generality, we have assumed g(t) = 0 on the inflow boundary j=0. 
This eliminates the influence of j=0 in the boundary polynomial, and reduces the order of 
the polynomials by one. Solving the resolvent equation for S and substituting into eqns (54, 
55, 56, 57) yields polynomials in k of the form 


( k 4 — 5 k 3 + 10k 2 — 9k + 9) 
(k 2 + 4k + 1) 

(4 tfc ) 

(58) 

(k 3 — 4k 2 + 5k — 8) 
(k 2 + 4k + 1) 

(3 rd ) 

(59) 

(k 2 - 2k + 7) 
(k 2 + 4k + 1) 

(2 nd ) 

(60) 

(k 2 — 2k — 11) 
(k 2 + 4k + 1) 

(l s< ) 

(61) 



(62) 

the roots of which are 



k = 2.286 ± 1.215», 0.2134 ± l-lS8i 

: (4 lfc ) 

(63) 

k = 3.218,0.3906 ± 1.527* 

(3 rd ) 

(64) 

k = 1 ± y/fti 

(2 nd ) 

(65) 

k = 1 ± 2y/3i 

(1 st ) 

(66) 


where i = y/~ T. It is apparent that |k| > 1 in all of these expressions. There are no 
eigenvalues or generalized eigenvalues, and the inflow boundary is stable to these closures. 
Given that the Cauchy problem and the two quarter-plane problems are stable, implies that 
the finite domain problem defined in eqn (12) is G-K-S stable for all boundary conditions 
specified thus far. 

A more rigorous example of the ability of the G-K-S theory to predict the stability of 
the fourth-order compact scheme is demonstrated by a pathological inflow boundary scheme. 
Since the inflow problem involves a NBS at the grid-point j=l which is biased in the “down- 
wind” direction, we suspect that it will be more sensitive to instability than the outflow 
boundary. We formulate a boundary scheme, which is a linear combination of the first- 
and second-order schemes, noting that they both have been shown to be stable boundary 
treatments. The resulting scheme which we shall use at the inflow boundary is 




= 2(1 + /?) 


— Vo + Vi 


where /3 is the amount of first-order influence in the formula. For /? = 0, the standard second- 
order closure is obtained. For all other values of /?, the scheme is formally first-order accurate. 
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We shall denote this scheme as (l J -4-4) since the inflow boundary is a one-parameter family 
of first order schemes, and the outflow is closed with a fourth-order formula. As was done 
previously, the dependence of the formula on the space derivative at j=0 is eliminated from 
eqn (67) by the inner-scheme written at j=l. The resulting expression is 

r« n , co\ (l + 2/3)<9V r 2 -(1 + 4/?)V 0 - 2(1 + fi)V\ + 3(1 + 2{3)V 2 , . 

I 2 (i + 2Z 3 ) - + 5 &T = 2Ai ' (68) 

Substituting Vj(t) = exp St <j>j, <f>j = 4> 0 k 3 into eqn (68) and noting the expression must be 
valid for the case g(t) = 0, results in a boundary scheme of the form 

([2(1 + 20) - i] + = -(1 + 0) + |(1 + 20)k. (69) 

Solving for S from the inner scheme eqn (42) and substituting it into eqn (69) produces a 
polynomial in n of the form 

(-1 + 2 /3)k 2 + (2 - 4/?)k - (7 + 22/3) = 0. (70) 


Solving for k yields 


k = 1 ± >/ 6 , 


(4/3 + 1) 
N (2/5- 1)’ 



(71) 


A double root exists for P = =£ and k = 1. Substituting k - 1 into the resolvent expres- 
sion yields 5 = 0, and a perturbation test shows that the boundary exhibits a generalized 
eigenvalue instability. Further inspection of eqn (71) shows that |/c| < 1 over the range 
X - P - = t' All other values of /3 need not be considered as candidates for instability since 
|/c| > 1. Substitution of the expression for k obtained from eqn (71), into eqn (69), yields 
an expression for S in terms of the parameter /?. Numerical evaluation of this expression 
shows that 1Z(S) > 0 for -.37 < (3 < =±. Thus, an eigensolution exists for this range, and 
the coupled inner-boundary scheme is unstable. 

To verify these findings, the model scalar wave equation described by eqns (28, 29, 30), 
was solved with the pathological inflow boundary conditions describe in eqn (67). Fourth- 
order boundary conditions were used at the outflow boundary. Figure (4) shows the results of 
the numerical investigation. Plotted is the logxo of the Li error of the solution integrated to 
a fixed time “T,” as a function of the parameter /3, ranging from [^-, |]. Two grid densities 
are shown in the study, and behave similarly. The theoretically predicted range of instability 
_.37 < ^ < zl i s replicated in the numerical study to within graphical limitations of the 
plot. In the unstable regime, the error grew exponentially with the number of iterations 
required to reach the time level T, and quickly became very large. 
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Although these boundary conditions are pathological, this study points to the fact that 
imposition of a first-order inflow boundary conditions is not sufficient to guarantee stability 
with the fourth-order compact inner-scheme. A similar experiment was performed on the 
outflow boundary using a linear combination of first- and second-order boundary conditions. 
For those boundary conditions, no eigensolutions could be found. That is not to imply that 
a completely arbitrary outflow boundary condition of at least first-order accuracy is always 
stable. It does indicate that the outflow is less susceptible to instabilities than the inflow 
boundary. 

Another note is appropriate, concerning the complexity of the G-K-S analysis for compact 
schemes. For the fourth-order compact inner-scheme, the polynomial of the highest degree 
which could not be factored (resulting in a numerical solution) was of degree four, and 
resulted from the closure with the fourth-order boundary conditions at the inflow plane. The 
G-K-S analysis of the fourth-order explicit scheme described by eqn (9) was performed by 
Strikwerda [8]. Stability polynomials of order eight were obtained for fourth-order closure 
at the inflow plane. It is apparent that the stability polynomials resulting from compact 
stencils are simpler expressions. This simplicity (and MACSYMA) will allow us to analyze 
schemes of sixth-order spatial accuracy without insurmountable algebraic polynomials at 
each boundary. 

6. SEMI-DISCRETE EIGENVALUE ANALYSIS 

The uniformly fourth-order explicit scheme (4, 4-4-4, 4) analyzed by Strikwerda [8] and the 
uniformly fourth-order compact scheme (4-4-4) presented here are both G-K-S stable for the 
semi-discrete problem and, therefore, will exhibit generalized stability for the fully discrete 
problem if advanced with a locally stable temporal scheme. This definition of stability ensures 
that the error of the numerical solution will remain uniformly bounded for all times by an 
exponentially growing amount. The exponential growth rate of the error is asymptotically (N 
— ► oo, where N is the total number of grid-points used) independent of the grid used. Thus, 
grid refinement studies with these methods, performed by integrating the governing equation 
to a fixed time level T on successively finer grids will demonstrate that the numerical solution 
converges to the exact solution at a rate of at least the order of the method. 

A disturbing feature of this stability definition is that the solution is not required to 
remain bounded for all times, even though the physical solution remains bounded for all 
times. Figure (5) shows a grid refinement study performed with the fourth-order compact 
(4-4-4) scheme demonstrating this behavior. The model equation was that described by 
eqns (28, 29, 30); the time interval was 0 < t < 100 and the grids used were 21, 41, and 81 
grid-points, respectively. Time was advanced with a fourth-order Runge-Kutta scheme in all 
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cases. The exact solution is a traveling sine wave of amplitude 1, for all times. Shown is the 
log w of the L 2 error, plotted as a function of time. Simulations on all three grids were run 
at various CFL’s. The initial portion of the simulation is characterized by nearly constant 
levels of error on all three grids. After a sufficiently long time, the unstable modes in the 
numerical solution dominate the spatial truncation error. From that point on, the solution 
diverges exponentially from the exact solution. Note that the growth rate in time, of the 
unstable modes of the solution is nearly grid independent, and that at any time T refining 
the grid by a factor of two results in a factor of 16 decrease in the error. It is also evident 
that at large times the actual error will be exponentially large. An interesting feature of 
the numerical method is that the exponential growth of the solution is dependent on the 
CFL used to advance the solution. For CFL = 1, the solution does not grow in time, while 
for CFL < a (a is some critical value less than CFL max ), exponential growth is observed. 
(This feature will be explained later, in terms of the amplification factor of the scheme.) 
Regardless of the CFL, fourth-order convergence is observed with the scheme. 

To understand the fundamental nature of the fixed grid, temporal divergence of the 
solution in the previous example, it is instructive to study the eigenvalue spectrum of the 
spatial discretization operator. As a semi-discretization, eqn (28) can be written in the form 
of eqn (5) as 

^ = MVj + j = 1, ...,N (72) 

where M is the N x N matrix describing the spatial discretization operator, and Bj g(t) 
represents the physical boundary data. Assuming that 

P~ 1 MP = S; P~ 1 U = V-, P~ 1 Bg(t) = H (73) 

where S is a diagonal matrix, and P~ l and P are similarity transforms composed of the left 
and right eigenvectors of the matrix M, respectively, eqn (72) takes the form 

^ = SU : + Hj, j = 1, (74) 

The solution to eqn (74) is 

Uj(t) = exp 5;< 17,(0) + f exp Sj 0 -T ) Hj(r)dT, j = 1, ... ,N. (75) 

J o 

In this form, the solution to eqn (74) depends exponentially on the eigenvalues Sj of the 
matrix M. Note that this solution assumes that the eigenvalues Sj are not degenerate, and 
that H(t) is not at a resonance frequency. If either of these situations occur, then the solution 
would include terms proportional to t p exy s P, where “p” is the order of the degeneracy. The 
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precise behavior depends on the temporal nature of H(t), but in general, for boundary data 
which remains bounded for all time, the solution grows only for the modes which have 
eigenvalues Sj with positive real parts. In addition, the growth rate will be governed by 
the eigenvalue with the largest positive component. Thus, any spatial discretization to the 
semi-discrete problem defined in eqn (72) will exhibit exponential divergence of the solution 
from the bounded physical solution, if it has an eigenvalue in the right half of the complex 
plane. 

Figures (6-7) show the eigenvalue spectrum resulting from the explicit fourth-order and 
the compact fourth-order spatial operators, closed at the boundaries with schemes of third- 
or fourth-order accuracy. In shorthand nomenclature, the explicit cases (3, 3-4-3, 3) and (4,4- 
4-4,4) are shown in Figure (6), and the compact cases (3-4-3) and (4-4-4) are shown in 
Figure (7). The spectrums are shown on grids of 21, 41, and 81 points, respectively. Note 
that closing the inner-schemes with third order NBS in both cases results in an eigenvalue 
spectrum which is bounded to the Left Half-Plane (LH-P), and that the uniformly fourth- 
order schemes cross over the imaginary axis into the Right Half-Plane (RH-P) of the complex 
plane. 

For long times, the maximum eigenvalues from the uniformly fourth-order schemes very 
accurately predict the exponential growth of the solution. In Figure (5), the solutions ob- 
tained from the (4-4-4) compact scheme grow exponentially in time. Assuming the er- 
ror can be represented functionally as e N (t) = e N (0)exp aN \ where N is the number of 
grid-points used in the spatial discretization, a growth rate am can be determined numer- 
ically. Similarly, from an eigenvalue determination, an effective grow rate as defined by 
ex P s |f?mox(At)| , can be calculated, where G max is the numerical amplification ob- 

tained from the temporal advancement scheme. For the fourth-order Runge-Kutta scheme 


G j — 1 T Sj At 4“ 


(Sj A* 2 ) (^Af 3 ) (S.At 4 ) 

2! 3! + 4! ’ 


7 = 1 , ...,N 


(76) 


and \G max \ will frequently correspond to the maximum eigenvalue 7 Z(S max ). Table (1) shows 
a comparison of the observed growth rate of the (4-4-4) compact scheme with that predicted 
from an eigenvalue determination. In each case, the maximum eigenvalue is used to predict 
the temporal grow of the solution. 


Grid 

® Numerical 

a K(Sma X ) 

21 

0.1321 

0.1315 

41 

0.1476 

0.1474 

81 

0.1537 

0.1479 


Table 1: Numerical vs. Theoretical Growth Rate; (4-4-4). 
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The agreement is very good, with a slight discrepancy in the comparison on the 81 grid-point 
case. In Figure (5), note the oscillatory nature of the growth of the solution. The uncertainty 
of the phase of the solution accounts for the discrepancy in the predicted growth rate in that 
case. 

A necessary condition for Lax stability of the finite domain semi-discretization can be 
expressed in terms of the eigenvalues of the spatial matrix operator as 

'R-(Sj) < w; u> > 0; j = 1, ••• , N (77) 

where Sj are the eigenvalues of the spatial operator, and N is an arbitrary number. It is 
apparent that the eigenvalue structure asymptotically approaches a bound in the RH-P as 

jV ► oo. All the fourth-order schemes presented thus far have satisfied this constraint. For 

the third-order NBS’s (3, 3-4-3, 3) and (3-4-3), the constant is a = 0, while for the fourth-order 
NBS’s, the constant is greater than zero. 

As mentioned earlier, a curious feature of the (4-4-4) (as well as other high-order spatial 
schemes) is that the growth of the solution is CFL dependent on all grids. For CFL’s close to 
the C F Lmaxi as determined from Von Neumann stability analysis, the schemes are bounded 
in time. For sufficiently small CFL’s, the schemes begin to diverge exponentially in time. 
(Again, it should be emphasized that in either case, the scheme is still G-K-S or Lax stable.) 
This behavior can be explained by noting a particular feature of the fourth-order Runge- 
Kutta time advancement scheme as well as some of the other locally stable time schemes. 
The stability bound of a time advancement scheme is defined as the locus of points in the 
complex plane where |z| < 1. Clearly, \z\ = 1 divides the plane into two regions. When 
the spatial eigenvalues (scaled by At) of a particular discretization lie entirely within the 
\Z\ = 1 boundary, the combined time-space scheme is generally stable. The stability regime 
of these schemes includes a semi-circular portion of the complex plane, centered at the origin 
and symmetric about the real axis and extending into the left half-plane. In addition, they 
contain a small part of the right half-plane, although not near the origin. If the spatial 
eigenvalues which lie in the RH-P are encompassed by the \Z\ = 1 boundary, the resulting 
scheme is stable. If a At is chosen such that the \Z\ = 1 line does not contain the RH-P 
eigenvalues, then the solution diverges with time. Figure (8) shows this feature for the (4-4- 
4) spatial scheme and the fourth-order Runge-Kutta scheme in time. For CFL’s which are 
near the CFL m ax , the maximum amplification rate \G ma x\ is less than one. For sufficiently 
small CFL’s, the \G max \ > 1 by an amount that is proportional to 1Z{S ma x), and the solution 
will diverge exponentially in time. It can be shown that a spatial scheme that has RH-P 
eigenvalues can always be made to diverge exponentially for sufficiently small CFL’s if a 
conventional third- or fourth-order R-K time advancement scheme is used. 
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7. ASYMPTOTIC STABILITY 

The previous discussion of Figure (5-8) brings out a subtle point in the Lax stability 
theory. In dealing with the numerical integration of time-dependent partial differential equa- 
tions, there are different limit processes that can be considered. One limit is the behavior of 
the numerical solution as the mesh size Ax -» 0 for a fixed time T. Another is the behavior 
of the solution for a fixed mesh as the time T tends to infinity. 

Stability, in the sense of Lax, addresses the first issue: boundedness of the numerical 
solutions as the mesh is refined at a fixed physical time. The essence of the Lax equivalence 
theorem is that if the numerical solution is bounded in this sense, then it converges to the 
true solution in the limit Ax — ► 0. To obtain an approximation to the true solution at time 
T, one integrates the IBVP up to time T on a sequence of grids as Ax -» 0. This sequence 
converges to the exact solution for all time levels T. 

Nothing in this definition excludes growth in time , and specifically allows exponential 
growth in time (see eqn (18)). Moreover, even if each of the quarter-plane problems is stable 
and allow no growth in time, the combined finite interval problem still allows exponential 
growth in time. (The Laplace transforms used in the G-K-S theory are legitimate only if 
growth in time is allowed.) 

Unfortunately, for genuinely time dependent problems, this stability definition might be 
too weak, in particular if long time integration is being carried out. The reason is that in 
order to achieve any reasonable accuracy for large times, one needs an excessive number of 
grid points. For long time numerical simulations to be useful, the solution of the semi-discrete 
problem defined in eqn (5), must be bounded in time as well. This means that for a fixed 
mesh N, the eigenvalues of the matrix M in eqn (5) have non-positive real part, and those with 
zero real part have (geometrical) multiplicity of 1. This is called Asymptotic Stability. 

An irksome feature of asymptotic stability is that, by itself, asymptotic stability does not 
imply Lax stability. There are numerous examples in the literature of fully discrete schemes 
that are asymptotically stable, but not Lax stable. The classic example of this is the case of a 
first order upwind spatial operator being advanced with an Euler explicit time advancement 
scheme. The eigenvalues for the fully discrete system are (1 - A), occurring a degenerate N 
times. Eigenvalue determination suggests that the CFL of the scheme should be 2, whereas, 
Von Neumann analysis and practical experience suggests that a CFL of one is the maximum 
stable CFL. The discrepancy can be explained by the fact that for 1 < CFL < 2, the matrix 
norm first grows rapidly before decaying asymptotically to zero. One might have suspected 
this difficulty by noting that in the semi-discrete case, the degenerate eigenvalues give rise to 
geometric growth in time, only later to be dominated by the exponentially decaying terms 
in the expressions. 
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An example for the semi-discrete case is now presented. Consider eqn (5) with g(t ) = 0. 
Let the matrix M be defined by the second-order central difference operator for the interior 
points. In matrix form, M can be written as 

'-11 

-10 1 0 

M = - ... (78) 

2 

0 -10 1 

-i 1 

where the boundary closures at grid-points j=l and j=N are artificially chosen. Using semi- 
discrete modal analysis, we determine a solution to eqn (5) by assuming the form V : (t) = 
exp st <t>j, with <f> = Ak 3 + B(=±y . The resolvent condition from the interior scheme yields 
Using the boundary conditions at point j = l and j=N to determine the values 
of A and B yields the expression for k of the form 

(* + l)((- 1 ) jv -k 2N ) = 0. (79) 

The roots to eqn (79) are k = -1 and k = iexp^ for j = l,N-l. The eigenvalues are thus: 

S = 0 and S = i cos for j=l,N-l, and are purely imaginary. The spatial discretization 
satisfies our definition of asymptotic stability for values of N, which are odd. 

The spatial discretization defined by eqn (78) admits a generalized eigenvalue instability 
at the inflow boundary. Using G-K-S analysis for the inflow boundary produces compatibility 
equations of the form 

2S k = k — 1 , j = 1 

2 S k 2 = k 3 - k, j = 2 (80) 

for which the only solution is k = 1, S = 0. The boundary condition is unstable to pertur- 
bations away from the unit circle and, therefore, exhibits a generalized eigenvalue instability 
at the boundary. It is apparent that asymptotic stability for the semi-discrete problem does 

not guarantee Lax stability. 

In addition to the previously mentioned examples, Reddy and Trefethen [14] have shown 
that it is not sufficient to consider the exact eigenvalues in determining the stability of a 
method. The famous Kreiss matrix theorem gives necessary and sufficient conditions for 
Lax stability in terms of the eigenvalues of the matrix M. A useful and equivalent test for 
determining stability is the analysis of the resolvent condition, which Reddy and Trefethen 
interpreted as involving not only the eigenvalues of the matrix M but also the e-pseudo- 
spectrum of the discretization matrix. This pseudo-spectrum is obtained by perturbing 
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the matrix M by an arbitrary matrix of norm e. Examples where the e-pseudo-eigenvalues 
determine the stability bounds of numerical methods are given by Reddy et al. [14], It is 
apparent that basing a stability definition on the eigenvalue determination of the spatial 
operator is not sufficient for stability. It is, however, useful for the present applications of 
higher-order scheme to restrict the allowable numerical discretizations to those possessing 
Lax stability and the property of bounded LH-P eigenvalues. For a broad class of spatial 
discretizations, these constraints are sufficient for the stability of the resulting numerical 
scheme. It is concluded that constraining the spatial operators to those possessing these 
properties is a useful and fairly general enhancement, and will be pursued in the remainder 
of this work. 

Before leaving fourth-order spatial discretizations, it is desirable to devise a uniformly 
fourth-order scheme. We already know that conventional discretization formulas at the 
boundaries result in G-K-S stable, but not asymptotically stable schemes for both the explicit 
and compact cases. The NBS’s used in each case relied on optimal order schemes at the 
boundary, where N+l constraints were used to devise the N th order scheme. Specifically, 
five grid-points in the explicit case, and four grid-points and one derivative condition in the 
compact case. If one removes the constraint of using optimal schemes at the boundaries, 
NBS with different dissipative characteristics can be devised and an asymptotically stable 
spatial scheme which is uniformly fourth-order can be found. 

We begin by devising an asymptotically stable fourth-order compact scheme, which we 
shall denote as (4 3 -4-4 3 ). The “4 3 ” signifies that the boundary point is closed with a fourth- 
order three-parameter family of schemes. The scheme defined at grid-point j=0 can be 
written as 

= (C 0 Uo + CiUi + C 2 U 2 + C 3 U 3 + C4U4 + C 5 U 5 + C 6 Ue + C 7 U 7 )/( Ax). (81) 

To be formally fourth-order accurate, Taylor series truncation analysis relates the coefficients 
in the following manner 

Go = —(a — 28/9 + 3227 + 13068)/5040 (82) 

Ci = +(a — 21/3 + 2957 + 5040)/720 

C 2 = -(a -26/3 + 2707 + 2520)/240 

C 3 = +(a — 25/9 + 2477 + 1680) / 144 

C 4 = —(a — 24/? + 2267 + 1260) /144 

C 5 = +(a- 23£ + 2077 + 1008)/240 

C 6 = —(a — 22/? + 1907 + 840)/720 

G 7 = +(a — 21/? + 1757 + 720)/5040 
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with similar expressions defined for the closure at the other end of the domain. By sys- 
tematically searching the three-parameter space spanned by the parameters a,/?, and 7, 
coefficients can be obtained which yield an eigenvalue spectrum which is bounded to the 
LH-P. The values of a,/?, and 7 are not unique, and no attempt to optimize the spectrum 
was made. A particular set of coefficients which make the scheme asymptotically stable are 
a = — 1560, /? = —355 and 7 = —35. Figure (9) shows the resulting spectrum from the 
(2-4-2), (3-4-3) and (4 3 -4-4 3 ) schemes. In all cases the eigenvalues are bounded to the LH-P, 
and the resulting scheme is asymptotically stable. 

Since the asymptotic stability condition 7Z(Sj) < u> for to = 0 is a very strong necessary 
condition, but not a sufficient condition, we still must show G-K-S stability for this case. 
For the outflow problem the model equation is 


dU_ 

dt 


Or r 

-7— = 0, x > 0, t > 0; 
Ox 


(83) 


no boundary condition is required in this problem, although a NBS is imposed at x = 0. 
Assuming a solution of the form Vj{t) = exp 5 * <^7, where <f>j = <f> 0 K 3 , and substitution into 
the inner scheme produces the resolvent condition for the eigenvalue S 

( — h 4 + k)S — 3( h k) (84) 

K K 

at each grid-point j > 1. At grid-point j=0, the parameter scheme produces an equation of 
the form 


3605 = -(+727 - 1370k + 1110k 2 - 875k 3 + 775k 4 - 552k 5 + 220k 6 - 35k 7 ). (85) 

Solving the resolvent condition for S and substituting into the boundary scheme yields a 
polynomial in k of the form 

(k - 1) 5 (35k 4 + 95k 3 - 168k 2 - 227k - 353) = 0. (86) 


Solving the polynomial for k produces no roots which are in magnitude less than one: thus, 
a stable condition. The possibility of k = 1 exists and must be checked for generalized 
eigenvalues. The condition is the same one tested for outflow stability in eqn (51), and was 
shown to be stable. Thus, the parametric fourth-order outflow scheme is G-K-S stable for 
the parameters a, f3 and 7 presented above. 

To show stability of the parameter scheme at the inflow, we study the partial differential 
equation 


dU 

dt 


r\ t t 

+ — — 0, x > 0, t > 0; 

ox 


(87) 


26 



with the boundary condition imposed at x = 0 of the form 


£/(0, t)=g(t), t> 0. 


(88) 


Eliminating ^ between the boundary scheme at j=0 and the inner-scheme at j=l, yields a 
combination boundary scheme of the form 


1440 


dUi 

dx 


+ 360 


dU 2 

dx 


— (+353t/ 0 + 1370C/! - 2190f/ 2 + 875 U 3 - 77567 

(89) 


+552f/ 5 - 220U 6 + 35t/ 7 )/(Ax) 


Substituting Vj(<) = exp S( <j> 0 K J into the inner-scheme and the boundary scheme produces 
the resolvent condition eqn (42), and 


(1440k + 360k 2 )S 


(+353C/ 0 + 1370k 1 - 2190k 2 + 875k 3 - 775k 4 + 552k 5 
-220k 6 + 35k 7 ). 


(90) 


Solving eqn (42) for S and substituting into eqn (91), with the condition that U 0 = 0 yields 
a polynomial in k of the form 


k(— 2950 + 2210k 1 - 2195k 2 + 1615k 3 - 1673k 4 + 1213k 5 - 293k 6 - 80k 7 + 35k 8 ) = 0. (91) 

Solving this polynomial results in no roots for which |k| < 1. The only possibility for instabil- 
ity is the condition |k| = 1. Again, this condition was checked previously for the other fourth- 
order schemes, and was shown to be stable for the inflow. The compact three-parameter 
family (4 3 -4-4 3 ) is stable for inflow and outflow, if the parameters are a = — 1560, /5 = —355 
and 7 = —35. In addition, the scheme produces an eigenvalue spectrum for the scalar wave 
equation which is bounded to the Left Half-Plane. It can be shown that treating the inflow 
and outflow boundaries of the explicit fourth-order scheme in a similar manner (4 3 , 4-4-4, 4 3 ) 
also produces an asymptotically stable scheme. Whether the resulting scheme is G-K-S 
stable was not checked in this work. 


8. SIXTH-ORDER SCHEMES 


As a last step in this work, the ideas and techniques used to analyze the fourth-order 
compact schemes are applied to sixth-order compact schemes. All the schemes tested here 
will be based on the sixth-order compact inner-scheme developed by Lele [15]. The scheme 
can be written in the form 


dUs - 1 + 3 du, 


dx 


dx 


+ 


dUj+\ _ —Uj-2 — 2817,-1 + 28{/,-+i + Uj + 2 


dx 


12Ax 


j =2, ... ,7V - 2. (92) 
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The scheme utilizes information from five-point explicitly and three point implicitly. As a 
consequence of the five point width of the stencil, NBS’s must be provided at two points at 
each end of the domain (j=0,l and j=N-l,N). The physical boundary condition is used at 
one of the inflow points. To ensure formal sixth-order accuracy for the hyperbolic problem, 
the boundary points must be closed with at least fifth-order formulas, which for the optimal 
schemes the short hand nomenclature would be (5, 5-6-5, 5). (Again, note that there is no 
ambiguity in comparing this nomenclature with that from the sixth-order explicit scheme, 
since it would require three boundary formulas at each end of the domain.) In keeping with 
the convention of this work, the closure at each end of the domain is done in a asymmetric 
manner, so that either the inflow or the outflow problem can be easily accommodated. 

It has proven extremely difficult to find G-K-S stable schemes which are formally sixth- 
order accurate. We, therefore, begin the discussion by presenting the stability analysis 
of a family of lower-order schemes. Two of the schemes in this family are the 1) (3,5-6- 
5,3) and 2) (4, 5-6-5, 4). The formal accuracies of these schemes are fourth- and fifth-order, 
respectively. The closure at the grid-point j=0 (with corresponding formulas written at j=N) 
is accomplished by 


dU Q 2 dUi 


dx 

dU o 


+ 3 


dx 

dUi 


— 5Uq + iUi -f- 1 / 2 
2Ax 

-nUo + 9U 1 +9U 2 -U 3 


(93) 

(94) 


dx dx 6Az 

for the third- and fourth-order schemes, respectively, while the fifth-order closure at the 
point next to the wall is accomplished in all cases by the scheme 


dUo dl\ „dU 2 -10Uo - 9Ux + 18U 2 + U 3 /0cO 

+ + = 3AJ • (95) 

Wider spatial stencils produce stability polynomials of dramatically increased complexity. 
Despite the fact that MACSYMA was used to determine all of the spatial formulas and 
the stability polynomial of the sixth-order schemes, the possibility for error still exists. 
Noting the ability of the G-K-S theory to accurately predict the stability envelop of the 
one-parameter family of fourth-order compact schemes (lM-4), a simple test was devised 
to verify the accuracy of the G-K-S calculations. A one-parameter family of schemes was 
created by combining the third- and fourth-order closure formula at each end of the domain. 
Symbolically, the combined scheme is represented by (3 1 , 5-6-5, 3 1 ), and is written as 


,au„ au, 

a|— 1- 2- 


dx 


dx 


— 5Uq T 4f/i + U 2 

2Ax 


] + (l 


dU 0 dUx 

«)[— + 3 


dx 


dx 


-ITUo + 9U X + 9 U 2 - Us 

6Ax 


(96) 
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For a — 0 (or 1), the scheme produces the optimal fourth-order (or third-order) variant, and 
for all other values of a, the formula is a third-order scheme. 

The model equation for the outflow quarter-plane problem is the same as described in 
eqn (83). Solutions of the form Uj(t) — exp St </> 0 satisfy the numerical scheme, giving the 
sixth-order inner-scheme the resolvent condition 

1 i 28 

( h 3 + k)S — ( — j 28/c + /c 2 )/ 12 (97) 

for the eigenvalue S. There are, in general, two roots to the resolvent equation which will 
yield |/c| < 1 for large TZ(S), since the resolvent is a polynomial of degree four. The other 
two roots will become exponentially unbounded as j becomes large and can be ignored. The 
general solution has the form Uj(t) = exp St (C X KF + C 2 n 2 j ). Substituting this expression into 
the boundary eqns (95,97) yields expressions for the constants C x and C 2 . The expressions 
are 


Ci-Fo(ki) + C 2 Fo(k 2 ) — 0 

CiFi(ki) + C 2 Fi(k 2 ) = 0 (98) 


where 

Fq(k) = (( — (6a + 18)/c — 18 /c 2 )5 — ((2a + 3) -(- (3a + 27)/c — (6a + 27)k 2 + (a — 3)/c 3 )) 
Fi(/c) = ((1 + 6« + 3« 2 )5 — (—10 — 9/c + 18 k 2 + 1k 3 )/3). (99) 


Equation (98) has only the trivial solution unless the determinant condition -Fo(ki)Fi(/C 2 ) — 
^i( k 2 )^o(ki) = 0 is satisfied. Solving the resolvent eqn (97) for S and substituting into the 
determinant condition yields an expression relating the two k's as 

(/ci — 1) («2 — 1) ((7a — 3)/ci k 2 + (—2a + 3)(«i + k 2 ) — 3a — 3) = 0. (100) 


Solving eqn (100) yields Kj = 1 or k 2 = 1, or 

(2a — 3 ) /c 2 T 3a -f- 3 
1 (7a — 3)« 2 — 2a + 3 

To obtain an additional independent relationship for /Cx and /c 2 , the resolvent condition eqn 
(97) is used. Solving the resolvent equation for S and noting that both k x and k 2 satisfy 
this expression for S, yields 





+ 28 /ci + /Ci^) 




28 

«2 


T 28/c 2 T k 2 2 ) 


1 2 (^ + 3 + K X ) 


+ 3 -f /c 2 ) 


( 102 ) 
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Combining eqn (101) and eqn (102) produces a single sixth-order polynomial in /c 2 , for which 
the roots can be found numerically. Equation (101) and eqn (97) are then used to determine 
the numerical values of and S , An eigensolution exists for the problem if, for |/ci| < 1 and 
|/c 2 | < 1, there exists a S with real part greater than zero. Solving the outflow polynomial 
for k 2 yields the result that the boundary is unstable for —9.16 < a < —1.86. The two 
limiting cases (a = 0 and ot = 1) for which the scheme is fourth- or third-order accurate, are 
stable. It should be noted that the possibility of instability also exists for the case /c a = /c 2 , 
and for K\ or k 2 = 1. None of these conditions showed instability, however. 

The model equation for the inflow quarter-plane problem is the same as described in eqns 
(87,88). Solutions of the form U 3 (t) = exp s< <f> 0 K j , will satisfy the numerical scheme, giving 
the sixth-order inner-scheme the resolvent condition 

(I + 3 + K )S = _(_1 - - + 2Sk + k 2 )/12. (103) 

k k k 2 k 

Eliminating ^ between the boundary schemes at grid-points j=0 and j=l, yields an ex- 
pression written for grid-point j = l of the form 

(6a +18)^- + 18^- = ((2a + 3)[/ 0 + (3a + 27)f/ 1 -(6a + 27)C/ 1 

(JX (JjC 

(104) 


+(o - 3)C 3 ))/(6 Ai). 


The general solution has the form Uj(t) = exp st (C\Ki J -\-C 2 K 2 i), Substituting this expression 
into the boundary expressions at j = 1 ,2 yields two expressions for the constants C\ and C 2 . 
The expressions are 

C\F\(k\) + C 2 F\(k 2 ) = 0 

C1F2M + C 2 F 2 (k 2 ) - 0 (105) 


where 

Fi(k) = (( — (6a + 18 )k — 18k 2 )5 — (+(3a + 27 )k — (6a + 27)/c 2 + (a — 3 )k 3 )) 

F 2 (k) = 1 . ( 106 ) 

The simple form of the expression F 2 = l results from reducing the modal equation at point 
j=2 (with Uq set to zero) by use of the resolvent condition (with U 0 not equal to zero). 
Equation (105) can only have a nontrivial solution if the determinant is identically zero. 
Solving eqn (103) for S and substituting into the determinant condition from eqn (105) 
yields 

(2a — 3)/ci 5 + (—5a + 15)/ci 4 — 30 ki 3 + (6a + 24)/ti 2 + (—22a — 33)/^! 1 — (a + 3) _ 

2«i 2 + 6K1 1 + 2 
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(2a - 3 )/c 2 5 + (—5a + 15 )k 2 4 - 30/c 2 3 + (6a + 24 )« 2 2 + (-22a - 33)^ - (a + 3) 

2k 2 2 + 6k 2 j + 2 


0 . 


(107) 

Together with eqn (102), they provide two equations for the unknown and n 2 . 

The inflow polynomial equations are far more difficult to solve since they do not factor 
appreciably. A change of variables from and k 2 to “x” and “y”, simplifies the algebra. 
Substituting 


«u + «2 = 2 y 
KiK 2 = x 


into eqn (102) yields 

(-8* - 8)?/ 2 + (— 12x 2 - 224x - 12 )y - 2x 3 - 166x 2 - 166x - 2 = 0 


(108) 


(109) 


which has the solution 

— (3x 3 + 56x + 3 ± V5x 4 + 2490x 2 + 5) , x 

y = - 4TT7 1 < u °) 

for x ^ —1. (The case x = -1 degenerates into y=0, for which /c x = — n 2 — ±1, a condition 
which produces no eigensolutions). Either root can be used since the final polynomial results 
from squaring an intermediate result to clear the square root in the expression. Substituting 
eqn (108) into eqn (107) and further simplifying with the expression for y from eqn (110) 
yields a ninth-order polynomial in the variable x, with coefficients that are functions of the 
variable a. Solving this expression numerically yields the roots for x. The values of y are 
then determined from eqn (110) and the values of /c x and n 2 are obtained from eqn (108). 
Again, the condition k = 1 and the case /c x = k 2 did not show instability over the parameter 
range tested in this study. For the inflow boundary condition, the instability envelop for the 
parameter a was determined to be -1.86 < a < -0.447. It is evident that the two degenerate 
conditions, a = 0 and a = 1, are G-K-S stable for the inflow. By combining the inflow and 
outflow results, the theoretical (G-K-S) range of instability for the one-parameter family of 
schemes (3 1 , 5-6-5, 3 1 ), was determined to be —9.18 < a < —0.447. The two degenerate cases, 
(3, 5-6-5, 3) and (4, 5-6-5, 4), were both G-K-S stable schemes. 

To determine the accuracy of the G-K-S calculations, another method of showing stability 
was used, and the results were compared. The necessary condition for stability on the 
eigenvalue structure ll(S m ax ) provides such a test. Figure (10) shows the results of an 
eigenvalue determination spanning the range -10<a<2, as determined numerically. The 
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maximum eigenvalue H(S max ) (the eigenvalue with the largest absolute real component), 
is plotted as a function of the parameter a for three grid densities. Since Lax stability in 
a finite domain requires that 7 Z(S) < u for w > 0, the 7l(S m ax ) should remain bounded 
with increasing grid density if the scheme is to be stable for a particular value of a. It is 
apparent from Figure (10), that the stability boundary at a = —9.15 is accurately predicted 
by G-K-S theory. The stability boundary at a = -0.45 is less well defined in the eigenvalue 
determination, and must be further investigated to show the correlation between G-K-S 
theory and eigenvalue determination. On the relatively coarse grids presented in Figure 
(10), the maximum eigenvalues near the limit a = —0.45 are growing with increasing grid 
density. Whether they are growing in a bounded manner determines if they satisfy the 
necessary condition for stability. Table (2) shows the behavior of TZ(S max ) for various grid 
densities at a = —0.40. 


Grid 

n(s max ) 

A% 

21 

1.000 

NA 

41 

1.250 

25 

81 

1.360 

8 

161 

1.582 

16 

321 

1.704 

8 

641 

1.780 

5 


Table 2: Grid convergence of 7 Z(S max ) for a — —0.40. 

As the number of grid-points becomes larger, the maximum eigenvalue asymptotes to a 
constant uq thus a stable condition. Note that the convergence to the asymptotic limit 
is slow for a = —0.40. Similar grid refinement studies at values of a = -0.45 and -0.50 
showed linear growth for all grids, thus an unstable condition. Based on a numerically 
determined eigenvalue determination over the range of —10 < a < 2, the G-K-S theory is 
accurately predicting the stability envelop. It should be noted that the slow convergence 
to the asymptotic limit (and the fact that the test is a necessary not sufficient one) makes 
testing for stability by numerical eigenvalue determination somewhat unreliable. For many 
non-borderline cases, it does provide an accurate measure of stability. 

The eigenvalue determination provides information on the asymptotic stability of the 
schemes as well. If for all the grids, the 7 £(S mar ) < 0, the scheme is asymptotically stable. 
For values of the parameter a < -9.15 and a > 0.4, eigenvalue determination indicates 
asymptotic as well as Lax stability of the resulting scheme. Figure (11) shows the eigenvalue 
spectrum of the (3, 5-6-5, 3) and (4, 5-6-5, 4) schemes. It is apparent that both satisfy the 
necessary condition of Lax stability, and that the (3, 5-6-5, 3) scheme is asymptotically stable. 
Because of the eigenvalues in the RH-P, the (4, 5-6-5, 4) scheme does not exhibit asymptotic 
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stability. Figure (12) shows a plot of the error of the solution to the scalar wave equation 
defined by eqns (28, 29, 30) when discretized with the (4, 5-6-5, 4) scheme. The time interval 
was 0 < t < 100 , and time was advanced with a fourth-order Runge-Kutta scheme. The 
logio of the Z /2 error is plotted as a function of time for three grid densities: 21, 41, and 81 
points, respectively. For a CFL = 1, the error does not grow in time. For CFL’s of order 
0.1, the error is nearly uniform for a period of time, then begins to grow exponentially with 
time. In all cases, doubling the grid decreases the error of the simulation by a factor of 16 
- 32 (error is dominated by either the fourth-order time, or the fifth-order space truncation 
terms.) The amplification is accurately predicted by the eigenvalue determination. Table (3) 
shows the numerical amplification rate compared with the Tt(S max ) , for which the agreement 
is excellent. 


Grid 

& Numerical 

a n(s maI ) 

21 

0.1245 

0.1228 

41 

0.1402 

0.1381 

81 

0.1351 

0.1354 


Table 3: Numerical vs. Theoretical Growth Rate; (4, 5-6-5, 4). 

The preceding examples of sixth-order schemes show the strength of G-K-S analysis to ac- 
curately predict the stability of complex higher-order schemes. They also show the intimate 
relationship between the eigenvalues of the spatial operator and the stability of the resulting 
scheme. 

We now present schemes that are formally sixth-order accurate, being closed at the 
boundaries by at least fifth-order stencils. Our first attempt is with optimal fifth-order 
closure at the boundaries, resulting in the scheme (5, 5-6-5, 5). Figure (13) shows the error 
of the simulation of the scalar wave equation defined by eqns (28, 29, 30). The behavior of 
this scheme is fundamentally different from the (4, 5-6-5, 4) scheme in several ways. On all 
grids, the error is always monotonically increasing in time. For CFL’s near the theoretical 
maximum value, the error increases at a lower rate, but is not suppressed as it was with 
the lower-order schemes. In addition, the exponential growth rate of the error increases 
with increasing grid density. On the grids shown, the error in the solution could not be 
systematically reduced by refining the grid and repeating the calculation to a specified time 
level T. In spite of these differences, the eigenvalue determination still accurately predicts 
the growth of the solution. Table (4) shows a comparison of the numerical and theoretical 
amplification rates. The theoretical values are determined from the 'JZ(Smar) for each grid. 
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Grid 

^ Numerical 

a n(s mai ) 

21 

0.7121 

0.7138 

41 

1.104 

1.010 

81 

1.749 

1.742 


Table 4: Numerical vs. Theoretical Growth Rate; (5, 5-6-5, 5). 

Again, the agreement is excellent. The solution grows at a rate which for long times is 
dominated by the eigenvalue with the maximum real part. The eigenvalue determination 
accurately predicts the behavior of the numerical solution even for this case which appears 
not to be Lax stable. 

Figure (14) shows a plot of the eigenvalue spectrum for the (5, 5-6-5, 5) scheme on a 
21, 41 and 81 grid. The TZ(S max ) is obviously increasing for these grids, and appears to 
be increasing without bound, as opposed to an asymptotic limit. This would violate a 
necessary condition for Lax stability. As was seen in the (3 1 , 5-6-5, 3 1 ) example, one cannot 
draw a precise conclusion from grid refined eigenvalue determination, although, the trends 
have in the cases presented thus far been the same. We must ultimately rely on G-K-S 
stability theory for the final answer as to whether the scheme is stable. 

We begin by determining the stability of the outflow boundary for the (5, 5-6-5, 5) scheme. 
The quarter-plane problem appropriate for this analysis is that described in eqn (83). No 
boundary conditions are necessary, but NRS’s are used at grid-points j=0,l. The closure at 
grid-point 0 is accomplished with 

dU 0 A dU x — 37C/ 0 + 8tA + 36*7 2 - SU 3 + U A , 11lX 

+ = iiAi <m) 

while that of grid-point j=l is accomplished with eqn (95). Solutions of the form Uj(t) = 
exp s< <f> oK,i, will satisfy the numerical scheme, giving the sixth-order inner-scheme the resol- 
vent condition shown in eqn (97). There will be two roots to the fourth-order polynomial in k 
which are |/c| < 1, and the general solution will have the form Uj(t) = exp^Ci/c^ + C 2 « 2 *). 
Substituting this expression into the two boundary conditions yields boundary expressions 
for the constants C\ and C 2 of the form 


C x F t 0 (ki) + C , 2 F o (« 2 ) = 0 

C 1 F 1 (k 1 ) + C 2 F 1 (k 2 ) = 0 (112) 

where 

F 0 (k) = ((1 +4ac x )5- (-37 + 8 * + 36/c 2 - 8 / c 3 + / c 4 )/ 12 ) 

Fi(k) = ((1 + 6/c + 3k 2 )S — ( — 10 — 9 k + 18« 2 + k 3 )/3). (113) 
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Equation (112) can only have a nontrivial solution if the determinant condition F 0 (ki)Fi(k 2 )~ 
F 1 (k 3 )F 0 (k 1 ) = 0 is satisfied. Solving the resolvent condition eqn(97) for S and substituting 
into the determinant condition produces 


- 0. 


(114) 


(115) 


(/Cl — 1) 4 (/C2 — 1) 4 («2 — «l)(4/C 2 Kl + K 2 T Ki — 6) _ 

144/ci/c 2 (/ci 2 “b 3/ci T 1)(^2^ "b 3 /c 2 T 1) 

Solving eqn (114) yields K\ = 1 or k 2 = 1, or Ki = k 2 , or the expression 

_ /c 2 - 6 

Kl Ak 2 + 1 

The first two roots are the same roots that have been shown previously to be stable. The 
condition that both roots must be |/c| < 1 precludes the last root. The third root can 
be shown to be stable by testing the derivative condition of the polynomial as follows. 
Multiplying and dividing eqn (112) by the non-zero rows and columns does not change the 
roots of the determinant conditions. The resulting expression is 

Tq(/C 2) — -fo(/Ci) 


To(/Ci) 


F i(*i) 


K 2 — Kl 

Fi(/c 2 ) - Fi(ki) 


= 0. 


K 2 - K 1 

Taking the limit as /ci —* k 2 yields the expression for the determinant condition F 0 (n)iS£i- 
F^) dF ^ = 0 . The resulting expression for /c is 

1)12 


(« 

144/c 2 (/c 2 + 3/c + l) 2 


0. 


(116) 


The outflow boundary is, thus, stable for the (5, 5-6-5, 5) scheme. 

The model equation for the inflow quarter-plane problem is the same as described in 
eqn (87, 88). Solutions of the form Uj(t) = exp 5 ‘^ 0 « j , will satisfy the numerical scheme, 
giving the sixth-order inner-scheme the same resolvent condition as was given in eqn (103). 
Again, the general solution will have the form Uj{t) = exp^Ci/Ci^ + C 2 k 2 j ). Substituting 
this expression into the two boundary conditions and making the simplification that Uq = 0 
gives two equations for the constants Ci and C 2 of the form 


CifiM T C 2 Fi(k 2 ) — 0 

CiT 2 (/ci) + C 2 F 2 (k 2 ) = 0 (117) 


where 

Fi(k) = ((2k + 3k 2 )S - (+44/c - 36/c 2 - 12/c 3 + /c 4 )/12 

F 2 (k) = 1. (118) 
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The determinant condition for which there is a nontrivial solution simplifies to Fi(ki) - 
Fi(k 2 ) — 0 and, along with eqn (102), provide two equations for the two unknowns sq and 
k 2 . Using the change of variables k 4 + « 2 = 2 y and k 4 k 2 = x, eqn (118) becomes 

—64 y 5 + (192 - 96 x)y 4 + (-16x 2 + 352x - 240)j/ 3 + (120x 2 - 504x + 160)y 2 

+ (8x 3 - 216x 2 + 360a: - 56 )y - 18x 3 + 142x 2 - 142x + 18 = 0. (119) 

Substituting the functional relationship y=y(x) provided by eqn (110) into eqn (119) and 
simplifying yields 

(x — l)(x + l) 5 (x 12 + 261X 11 + 24298x 10 + 864903x 9 + 864903x 9 + 555871 lx 8 
+ 16502410x 7 + 27479264x 6 + 28538822x 5 + 5255107x 4 + 429169x 3 

+ 18614x 2 + 435x + 5) = 0. (120) 

Two roots to this polynomial give rise to eigenvalues S which are in the RH-P. They are x = 
(-0.01969168445 ±0.02398022319*), which yields k x = (0.157055621 +0.943601129*), with 
k 2 = (—0.02810826491 , +0.01619023438*) and S = (0.0428389 ± 1.39944*). The numerical 
solutions satisfy the governing equations to approximately machine precision (1.0e-13). Thus, 
the inflow for the (5, 5-6-5, 5) scheme is G-K-S unstable. This verifies the trends presented 
earlier in the eigenvalue grid refinement analysis, and the simulation of the scalar wave 
equation, both of which indicated the sixth-order scheme was Lax unstable. 

The unstable (5, 5-6-5, 5) scheme previously discussed was implemented using optimal 
fifth-order boundary formulas. By relaxing the constraint of optimal order schemes at the 
boundaries, the possibility of a fifth-order closure which will be G-K-S stable still exists. 
Taylor series truncation analysis was used to develop parametric relations for the closure 
formulas at the two NBS’s. They were constrained to be fifth-order, and explicit in nature. 
To facilitate a wide range of closures, each point was given two degrees of freedom. The 
symbolic formula for the new scheme is, thus, (5 2 ,5 2 -6-5 2 ,5 2 ). The scheme defined at grid- 
point j=0 and 1 can be written as 

qrr 

= (Co^o + C\U\ + C2U2 + C 3 U 3 + C4U4 + C5O5 + CeUe + CjUr) / (Ax) (121) 

= ( D 0 U 0 + D 1 U 1 + D 2 U 2 + D 3 U 3 + D 4 U 4 + D 5 U 5 + D 6 U 6 + D 7 U 7 )/( Ax) (122) 

where 

Co = — (c*o — 28/3 0 + 13068)/ 5040 (123) 
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Cl = +(a 0 - 27ft + 5040)/720 
C 2 = — (a 0 — 26ft + 2520) /240 
C 3 = +(a 0 - 25ft + 1680)/144 
C 4 = — (a 0 — 24ft + 1260)/144 
C 5 = +(a 0 - 23ft + 1008)/240 
C e = ~(a 0 - 22ft + 840)/720 
C 7 = +(a 0 - 21ft + 720)/5040 

and 

Do = -(c*i - 21ft + 720)/5040 (124) 

Di = +(cti - 20ft - 1044)/720 

D 2 = -{on- 19ft - 720)/240 

£) 3 = +( ai - 18/?! - 360)/144 

Z) 4 = — (c*i - 17A - 240)/144 

7^5 = +(a! - 16A - 180)/240 

£>6 = -(«i - 15A - 144)/720 

D 7 = +(ai - 14ft - 120)/5040 

with similar expressions defined for the closure at the other end of the domain. By systemat- 
ically searching the four- parameter space spanned by the parameters a 0 ,/?o,ai, and ft with 
an eigenvalue code, an arbitrary eigenvalue spectrum can be approximated. A particular 
set of coefficients for which the eigenvalue spectrum is bounded to the Left Half-Plane is 
a 0 = 1809.257, ft = -65.1944,0! = -262.16 and ft = -26.6742. The values are not unique 
and no attempt has been made to find optimal values of these coefficients. The eigenvalue 
spectrum for this case is shown on Figure (14). Note that the shape of the spectrum is similar 
to that of the (5, 5-6-5, 5) scheme, but that the TZ(S max ) < 0 instead of increasing without 
bound. The scheme satisfies the necessary condition for Lax stability, and is asymptotically 
stable by our definitions. 

The stability analysis for the two quarter-plane problems involved in establishing the 
G-K-S stability of the new (5 2 ,5 2 -6-5 2 ,5 2 ) is extremely formidable. It pushes MACSYMA to 
the limits of its capabilities on present machines. In addition, 128-bit arithmetic is required 
to ensure precision when determining the roots of the resulting polynomials in “x”. The 
scheme has been shown to be G-K-S stable for the parameters given above for the inflow and 
outflow problems. Thus a formally sixth-order scheme has been developed which is G-K-S 
(Lax) stable, and asymptotically stable for the scalar case. 
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To verify its accuracy, Table (5) shows a grid refinement study performed with the new 
sixth-order scheme. The model problem is the scalar wave equation defined by eqns (28, 29, 
30). The time advancement scheme is the fourth-order Runge-Kutta algorithm, with a CFL 
of 0.10. Temporal refinement studies were performed to ensures that the leading error terms 
on all grids were from the spatial discretization operator. Listed is the grid and its error, 
and the slope between each successive refinement. The asymptotic stability of the spatial 
operator ensure that the solution does not grow exponentially for long times. 


Grid 

L 2 error 

slope 

21 

-2.363 

NA 

31 

-3.582 

-6.9 

41 

-4.225 

-5.1 

51 

-4.724 

-5.1 

61 

-5.150 

-5.4 

81 

-5.849 

-5.6 

101 

-6.406 

-5.7 

121 

-6.867 

-5.8 


Table 5: Grid refinement study of the (5 2 , 5 2 -6-5 2 , 5 2 ) scheme. 

The data point corresponding to 21 grid-points is erroneous because the grid is too coarse 
for the scheme to exhibit its higher-order properties. Note that in this example the scheme 
becomes at least fifth-order accurate at approximately 10 grid-points/ (27T radians). In the 
limit N — ► oo where N is the number of grid-points, the scheme shows a slope of -6, the 
formal accuracy of the inner-scheme. Note that for N small, the four points which are treated 
with fifth-order accuracy degrade the formal accuracy by one degree. 

An asymptotically stable sixth-order spatially accurate scheme has been developed for 
use in a method-of-lines discretization of a hyperbolic partial differential equation. The 
eigenvalues of the new scheme are for the scalar case bounded to the Left Half-Plane of the 
complex plane for all N. The necessary condition for Lax stability is, therefore, satisfied. 
In addition, the scheme has been shown to be G-K-S stable for the combined inflow and 
outflow quarter-plane analysis, and is, therefore, formally Lax stable for the scalar case. For 
the hyperbolic system of equations, the use of any of the Lax stable schemes presented in this 
work guarantees the Lax stability of the resulting spatial discretization if the boundaries are 
imposed in characteristic form. The concept of asymptotic stability does not carry over from 
the scalar case to the system. For the case of the system with all of the physical eigenvalues of 
the same sign, the asymptotic stability is still retained. For the case of mixed eigenvalues, it 
is easy to show that even though being asymptotically stable for the scalar case, exponential 
growth of the solution may occur for boundaries that are imposed in characteristic form. 
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Work continues in this area, to determine a stronger necessary condition for asymptotic 
stability that will allow the use of scalar analysis to determine spatial schemes which are 
asymptotically stable for the system. 

9. CONCLUSIONS 

The stability characteristics of various compact fourth- and sixth-order spatial operators 
were assessed using the theory of Gustafsson, Kreiss and Sundstrom (G-K-S) for the semi- 
discrete IBVP. The class of central difference interior schemes with asymmetrically closed 
boundaries was analyzed. Because of formal accuracy considerations, those schemes with 
boundary closures of at least ( N — l) </l spatial order for an (N) th order inner scheme were 
the focus of the work. It was found that conventional third- or fourth-order boundary 
conditions, when coupled with the fourth-order compact inner-scheme, resulted in a G-K : S 
stable scheme. For the sixth-order compact inner-scheme, the conventional boundary closures 
of fifth- and higher-order were found to be G-K-S unstable. Fourth-order and lower-order 
closure formulas were found to be G-K-S stable. These results were then generalized to the 
fully discrete case using a recently developed theory of Kreiss, which states that under weak 
constraints, the stability of the semi-discrete operator implies stability of the fully discrete 
operator if a locally stable temporal method is used. 

The conventional definition of stability was then sharpened to include only those spatial 
discretizations that are asymptotically stable (bounded Left Half-Plane eigenvalues). Many 
of the higher-order schemes which are G-K-S stable were found to not be asymptotically 
stable. Fourth-order boundary conditions were found to be asymptotically unstable for the 
schemes, tested, specifically: 1) (4-4-4), and 2) (4, 5-6-5, 4). A series of compact fourth- and 
sixth-order schemes which were both asymptotically and G-K-S stable were then developed. 
The constraint of optimal accuracy from a specific number of constraints was abandoned, 
thus enabling several- parameter family boundary closures to be developed. A three param- 
eter family uniformly fourth-order scheme (4 3 -4-4 3 ), as well as a four-parameter sixth-order 
scheme with fifth-order boundaries (5 2 , 5 2 -6-5 2 , 5 2 ) were developed which were asymptotically 
stable. No attempt was made to optimize the parameters. 

All the schemes which were found to be G-K-S stable were subjected to extensive com- 
parisons between the G-K-S stability predictions, semi- discrete eigenvalue determination and 
numerical simulations. In all cases, consistent and complementary results were achieved with 
all three methods. In addition, it was shown that the eigenvalue determination accurately 
predicted the exponential divergence of the solution for the cases which were not asymptot- 
ically stable. Work continues in developing asymptotic definitions for the case of systems of 
equations. 
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PADE BOUNDARY CONDITION ACCURACY 



Fig-1 Comparative study of the effects of boundary closure on the formal accuracy of the 
fourth-order compact inner scheme. 
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PADE BOUNDARY CONDITION ACCURACY 



Fig-2 Comparative study of the effects of outflow boundary closure on the formal accuracy 


of the fourth-order compact inner scheme. 
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PADE BOUNDARY CONDITION ACCURACY 



YO 1.2 i A ' .6 ' .8 

Log(Grid Points) 


Fig-3 Comparative study of the effects of inflow boundary closure on the formal accuracy 
of the fourth-order compact inner scheme. 
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0.0 


STABLE INFLOW CONDITIONS 



Fig-4 Numerical determination of the Z. 2 error at a time level “T” resulting from a one- 
parameter family of closure formulas at the inflow boundary. The inner scheme is the 
fourth-order Pade scheme. 
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ERROR OF D ADE SCHEME (4-4-4) 



Fig-5 The L 2 error as a function of time, resulting from the uniformly fourth-order Pade 
scheme (4-4-4), for various CFL’s and grid densities. 
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EIGENVALUE SPECTRUM (4-th LINEAR) 
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Fig-6 The numerically determined eigenvalue spectrum resulting from the explicit fourth- 
order inner scheme, closed with either third- or fourth-order boundary schemes. 
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EIGENVALUE SPECTRUM (4-th PADE) 
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Fig-7 The numerically determined eigenvalue spectrum resulting from the fourth-order Pade 
inner scheme, closed with either third- or fourth-order boundary schemes. 
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AMPLIFICATION FACTOR (4 th PADE ) 



CFL 


Fig-8 The maximum numerical amplification factor as a function of CFL, plotted for the 
uniformly fourth-order Pade scheme, with a fourth-order Runge-Kutta time advance- 
ment scheme. 
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EIGENVALUE SPECTRUM (NTH PADE) 
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Fig-9 The numerically determined eigenvalue spectrum resulting from the fourth-order Pade 
inner scheme, closed with either second-, third- or newly developed fourth-order bound- 
ary schemes, on various grids. 
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MAXIMUM EIGENVALUE (p,5-6-5,p) 



alpha 


Fig-10 The maximum real part of the numerically determined eigenvalue spectrum, result- 
ing from a sixth-order Pade inner scheme and a one-parameter family of third-order 
boundary schemes. The ordinant a is the boundary parameter. 
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EIGENVALUE SPECTRUM (6 th PADE) 


a (3,5-6-5,3)21 grdpts 
O (3,5-6-5,3)41 grdpts 
o (3,5-6-5,3)81 grdpts 
+ (^,5-6-5,4)21 grdpts 
* (^,5 — 6 — 5,4)4 1 grdpts 
x (^,5-6-5, 4)8 ' grdpts 


Fig-11 The numerically determined eigenvalue spectrum resulting from the sixth-order Pade 
inner scheme, closed with either third- or fourth-order boundary schemes, on various 





ERROR OF PACE SCHEME (4,5-6-5.4) 



Fig-12 The L 2 error as a function of time, resulting from the uniformly sixth-order Pade 
scheme with fourth-order boundary closure (4, 5-6-5, 4), for various CFL’s and grid 
densities. 
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ERROR OF PADE SCHEME (5, 5-6-5, 5) 



Fig-13 The T 2 error as a function of time, resulting from the uniformly sixth-order Pade 
scheme with fifth-order boundary closure (5, 5-6-5, 5), for various CFL’s and grid den- 
sities. 


54 



EIGENVALUE SPECTRUM (6 th PADE) 


A (5,5 6 5,5)21 grdpts 
O (5,5-6-5,5)41 grdpts 
O (5,5 -6-5,5)81 grdpts 
+ (5 2 ,5 2 -6-5 2 ,5 2 )21 grdpts 
* (5 2 ,5 2 -6 — S 2 ^ 2 )'! 1 grdpts 
x (5 2 ,5 2 -6-5 2 ,5 2 )81 grdpts 


Fig-14 The numerically determined eigenvalue spectrum resulting from the sixth-order Pade 
inner scheme, closed with conventional fifth- or newly developed fifth-order boundary 
schemes, on various grids. 








16. Abstract 

The stability characteristics of various compact fourth- and sixth-order 
spatial operators are assessed using the theory of Gustafsson, Kreiss and Sundstrom 
(G-K-S) for the semi-discrete Initial Boundary Value Problem (IBVP) . These results 
are then generalized to the fully discrete case using a recently developed theory 
of Kreiss. In all cases, favorable comparisons are obtained between G-K-S theory, 
eigenvalue determination, and numerical simulation. The convential definition of 
stability is then sharpened to include only those spatial discretizations that are 
asymptotically stable (bounded, Left Half-Plane eigenvalues) • It is shown that many 
of the higher-order schemes which are G-K-S stable are not asymptotically stable. 

A series of compact fourth- and sixth-order schemes, which are both asymptotically 
and G-K-S stable for the scalar case, are then developed. 
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